#ifndef OPTION1D_HPP
#define OPTION1D_HPP

#include <iostream>
#include "Eigen/Dense"
using namespace std;
using namespace Eigen;

enum BoundCondiType
{
    Dirichlet,
    Neumann,
    Mixed
};

enum MultigridMethod
{
    Vcycle,
    FMG
};

namespace Operators
{
    enum Restriction
    {
        Injection,
        FullWeighting
    };

    enum Interpolation
    {
        linear,
        qudratic
    };
}

MatrixXd RestrictionMatrix(Operators::Restriction op, int n) // n = 2^i
{
    MatrixXd m = MatrixXd::Zero(n / 2 + 1, n + 1);
    if (op == Operators::Restriction::Injection)
    {
        for (int i = 0; i < n / 2 + 1; i++)
        {
            m(i, 2 * i) = 1;
        }
    }
    else if (op == Operators::Restriction::FullWeighting)
    {

        m(0, 0) = 1;
        m(n / 2, n) = 1;
        for (int i = 1; i <= n / 2 - 1; i++)
        {
            m(i, 2 * i) = 0.5;
            m(i, 2 * i - 1) = 0.25;
            m(i, 2 * i + 1) = 0.25;
        }
    }
    // cout << m << endl;
    return m;
}

MatrixXd InterpolationMatrix(Operators::Interpolation op, int n)
{
    MatrixXd M = MatrixXd::Zero(2 * n + 1, n + 1);
    if (op == Operators::Interpolation::linear)
    {
        for (int i = 0; i < n + 1; i++)
        {
            M(2 * i, i) = 1;
        }
        for (int i = 0; i < n; i++)
        {
            M(2 * i + 1, i) = 0.5;
            M(2 * i + 1, i + 1) = 0.5;
        }
    }
    else if (op == Operators::Interpolation::qudratic)
    {
        M(0, 0) = 1;
        for (int i = 1; i <= n - 1; i++)
        {
            M(2 * i, i) = 1;
            M(2 * i - 1, i - 1) = 0.375;
            M(2 * i - 1, i) = 0.75;
            M(2 * i - 1, i + 1) = -0.125;
        }
        M(2 * n - 1, n) = 0.375;
        M(2 * n - 1, n - 1) = 0.75;
        M(2 * n - 1, n - 2) = -0.375;
        M(2 * n, n) = 1;
    }
    //cout << M << endl;
    return M;
}

#endif